function test_functions()

  
  params = h2params(1, 10, 32, 1, 512*4*2, 1);
  load('nf.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');



  z = params.Cz';
  T = 1;

  
  % -- z^2
   % Fh2s6 = h2s6(z, 1);
   % Fh2s7 = h2s7(z, 1); 

   % h = heq(Fh2s6);
   % plot(z, h-hs6(z));

   % hold on

   % h = heq(Fh2s7);
   % plot(z, h-hs7(z));

   % hold off

   % plot(z, ks7(z))

   % hold off
   % k = keq(k2s3(z, 1));
   % plot(z, k-ks3(z));

   % hold on

   % k = keq(k2s5(z, 1));
   % plot(z, k-ks5(z));

   % k = keq(k2s6(z, 1));
   % plot(z, k-ks6(z));

   % k = keq(k2s7(z, 1));
   % plot(z, k-ks7(z));

   % hold off

   % j = jeq(j2v3(z, 1));
   % plot(z, j - jv3(z));
   % hold on
   % 
   % j = jeq(j2v4(z, 1));
   % plot(z, j - jv4(z));

   % j = jeq(j2v5(z, 1));
   % plot(z, j - jv5(z));

   % j = jeq(j2v7(z, 1));
   % plot(z, j - jv7(z));

  a = aeq(a2t2(z, 1));
  plot(z, a-at2(z), 'b')
  hold on

  a = aeq(a2t4(z, 1));
  plot(z, a-at4(z), 'r')

  a = aeq(a2t5(z, 1));
  plot(z, a-at5(z), 'g')

  hold off 
  % *********************
  %  diff eq functions
  % *********************

  function rheq = heq(h)
    % h -- z^2
    hz = chebgrad(h, params.Cz(1), params.Lz, 1);
    hz = 2.*z.*h + z.^2.*hz;
    hzz = chebgrad(hz, params.Cz(1), params.Lz);
    
    rheq = z.^3 .* (-2 .* hz + z.*hzz);
  end

  function rkeq = keq(k)
    % k -- 1/z^2 
    kz = chebgrad(k, params.Cz(1), params.Lz, 1);
    kz = -2.*k + z.*kz;  % 1/z^3

    rkeq = 2 .* k + 2.* kz;
  end

  function rjeq = jeq(j)
    % j -- 1/z
    jz = chebgrad(j, params.Cz(1), params.Lz, 1);
    jz = -j + z.*jz; %1/z^2

    jzz = chebgrad(jz, params.Cz(1), params.Lz, 1);
    jzz = -2.*jz + z.*jzz; %1/z^3

    rjeq = 0.5 .* z.^2 .*(4.*jz + jzz);
  end

  function raeq = aeq(a)
    az = chebgrad(a, params.Cz(1), params.Lz, 1);
    azz = chebgrad(az, params.Cz(1), params.Lz, 1);

    raeq = az./z + z.^2.*az./2 + 1./2.*(-1+z.^3).*azz;
  end





  % *************
  % source functions
  % *************

  function rhs6 = hs6(z)
    rhs6 = (-1/2).*z.^4;
  end

  function rhs7 = hs7(z)
    rhs7 = (1/3).*z.^4.*(1+z+z.^2).^(-2).*((-3)+(-6).*z+2.*3.^(1/2).*pi.*z+(-3).*z.^2+3.^(1/2).*pi.*z.^2+(-2).*3.^(1/2).*z.*(2+z).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*z.*(2+z).*log(1+z+z.^2));
  end

  function rks3 = ks3(z)
    rks3 = 2;
  end

  function rks5 = ks5(z)
    rks5 = 1/2;
  end

  function rks6 = ks6(z)
    rks6 = (-2)+(-1/2).*z.^3;
  end

  function rks7 = ks7(z) 
    rks7 = z.^3.*(1+z+z.^2).^(-1).*((-1)+4.*z.^(-3)+4.*z.^(-2)+2.*z.^(-1)+(1/3).*z.^(-4).*((-4)+(-4).*z+z.^3+z.^4).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
  end


  function rjv3 = jv3(z)
    rjv3 = (1/2).*z.^2;
  end

  function rjv4 = jv4(z)
    rjv4 = (-1/2).*z.^4.*(1+z+z.^2).^(-1);
  end

  function rjv5 = jv5(z)
    rjv5 = (1/2).*(1+z+z.^2).^(-2).*(z.^2+2.*z.^3+3.*z.^4+z.^5+(-1).*z.^6);
  end

  function rjv7 = jv7(z)
    rjv7 = (-1/2).*z.^4.*(2+z).*(1+z+z.^2).^(-2);
  end

  function rat2 = at2(z)
    rat2 = (1/3).*z.^(-1).*(3.^(1/2).*pi+(-6).*z.*(1+z).*(1+z+z.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2));
  end

  function rat4 = at4(z)
    rat4 = (-1/12).*z.^(-1).*(3.^(1/2).*pi+6.*z.^3.*(1+z+z.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2));

  end

  function rat5 = at5(z)
    rat5 = (1/2).*(((-2)+(-2).*z+z.^2).*(1+z+z.^2).^(-1)+(1/2).*z.^(-1).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
  end

%     for (iz = 1:params.Nz)
%       x = z(iz);
%       if (x < 0.001)
%         rks3(iz) = 
%       else
%         rks3(iz) = 
%       end
%     end



    % *********************
    %   h - z^2
    % *********************

    function resh2s6 = h2s6(z, T) % 1/18/13
      % has z.^2 factored out. but only z, no T (that goes for everybody)!!
      % -1/2r^4
      x = z(:).*T;
      resh2s6 = 0.*z(:) +  (T.^2)./4;
    end

    function resh2s7 = h2s7(z, T)
      % has z.^2 factored out. but only z, no T!!
      %resh2b = 0;
      x = z(:).*T;
      resh2s7 = chebinterp(nh2s7, x, 0, 2);
      resh2s7 = resh2s7 .* T .* T;
    end

    % *********************
    %   k - 1/z^2
    % *********************

    function resk2s3 = k2s3(z, T) % 1/18/13
      % k2s3 with 1/z^2 extracted
      % 2 S3
      resk2s3 = 0.*z(:) - 1;
      resk2s3 = resk2s3./(T.^2);
    end
    function resk2s4 = k2s4(z, T)
      % placeholder
      resk2s4 = 0.*z(:);
      resk2s4 = resk2s4./(T.^2);
    end
    function resk2s5 = k2s5(z, T) % 1/18/13
      % k2s5 with 1/z^2 extracted
      % 1/2 S5
      resk2s5 = 0.*z(:) - 1/4;
      resk2s5 = resk2s5./(T.^2);
    end
    function resk2s6 = k2s6(z, T) % 1/18/13
      % k2s6 with 1/z^2 extracted.
      % -(1+4r^3)/(2r^3) S6
      x = z(:).*T;
      resk2s6 = 1+(-1/8).*x.^3;
      resk2s6 = resk2s6./(T.^2);
    end
    function resk2s7 = k2s7(z, T)
      % k2s7 with 1/z^2 extracted.
      x = z(:).*T;
      resk2s7 = chebinterp(nk2s7, x, 0, 2);
      resk2s7 = resk2s7./(T.^2);
    end

    % *********************
    %   j - 1/z
    % *********************

    function resj2v3 = j2v3(z, T) % 1/18/13
      % j2v3 with 1/z extracted!
      % same as in mark's paper
      % 1/(2r^2) V3
      resj2v3 = 0.*z(:) + -1/2;
      resj2v3 = resj2v3./T;
    end

    function resj2v4 = j2v4(z, T) % 1/18/13
      % with 1/z extracted
      % same as in mark's paper!
      % -1/(2r^2(1+r+r^2)) V4
      resj2v4 = zeros(length(z), 1);
      for i=1:length(z)
        x = z(i).*T;
        if (x < 0.001)
          resj2v4(i) = (-1/4).*x.^2+(1/10).*x.^3+(-1/28).*x.^5+(1/40).*x.^6+(-1/70).*x.^8+(1/88).*x.^9+(-1/130).*x.^11+(1/154).*x.^12+(-1/208).*x.^14+(1/238).*x.^15+(-1/304).*x.^17+(1/340).*x.^18+(-1/418).*x.^20; % feb 27
        else
          resj2v4(i) = (1/54).*x.^(-2).*(9.*((-2)+x).*x+2.*3.^(1/2).*pi.*((-1)+x.^3)+(-12).*3.^(1/2).*((-1)+x.^3).*atan(3.^(-1/2).*(1+2.*x)));
        end
      end
      resj2v4 = resj2v4./T;
    end

    function resj2v5 = j2v5(z, T) % 1/18/13
      %j2v5 with 1/z extracted, the resulting thing has lowest order const.
      % (-1+r+3r^2+2r^3+r^4)/(2r^2(1+r+r^2)^2) (V5-V6)
      resj2v5 = zeros(length(z), 1);
      for i=1:length(z)
        x = z(i).*T;
        if (x < 0.001)
          resj2v5(i) = (-1/2)+(-1/10).*x.^3+(3/28).*x.^5+(-1/10).*x.^6+(3/35).*x.^8+(-7/88).*x.^9+(9/130).*x.^11+(-5/77).*x.^12+(3/52).*x.^14+(-13/238).*x.^15+(15/304).*x.^17+(-4/85).*x.^18+(9/209).*x.^20;
        else
          resj2v5(i) = real((1/54).*x.^(-2).*(9.*((-8)+x).*x+(-2).*3.^(1/2).*pi.*(4+11.*x.^3)+12.*3.^(1/2).*((4+3.*x.^3).*atan(3.^(-1/2).*(1+2.*x))+(sqrt(-1)*(-2)).*x.^3.*(log(sqrt(-1)+(-1).*3.^(1/2)+(sqrt(-1)*2).*x)+(-1).*log(sqrt(-1)+3.^(1/2)+(sqrt(-1)*2).*x)))) );
        end
      end
      resj2v5 = resj2v5./T;
    end

    function resj2v6 = j2v6(z, T) 
      resj2v6 = -1.*j2v5(z,T);
    end

    function resj2v7 = j2v7(z, T) % 1/18/13
      %j2c with 1/z extracted, the resulting thing has lowest order const.
      % (1+2r)/(2r(1+r+r^2)^2) V7
      resj2v7 = zeros(length(z), 1);
      for i=1:length(z)
        x = z(i).*T;
        if (x < 0.001)
          resj2v7(i) = (-1/2).*x.^2+(3/10).*x.^3+(-5/28).*x.^5+(3/20).*x.^6+(-4/35).*x.^8+(9/88).*x.^9+(-11/130).*x.^11+(6/77).*x.^12+(-7/104).*x.^14+(15/238).*x.^15+(-17/304).*x.^17+(9/170).*x.^18+(-10/209).*x.^20;
        else
          resj2v7(i) = (1/27).*x.^(-2).*((-9).*((-2)+x).*x+3.^(1/2).*pi.*(2+x.^3)+(-6).*3.^(1/2).*(2+x.^3).*atan(3.^(-1/2).*(1+2.*x)));
        end
      end
      resj2v7 = resj2v7./T;
    end




    % ***** % 
    %   a   %
    % ***** %

    % no extraction

    function resa2t2 = a2t2(z, T)
      x = z.*T;
      resa2t2 = chebinterp(na2t2, x, 0, 2);
    end

    function resa2t3 = a2t3(z, T)
      resa2t3 = a2t2(z,T);
    end

    function resa2t4 = a2t4(z, T)
      x = z.*T;
      resa2t4 = chebinterp(na2t4, x, 0, 2);
    end
    function resa2t5 = a2t5(z, T)
      x = z.*T;
      resa2t5 = chebinterp(na2t5, x, 0, 2);
    end

end
